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c/3 ' The compacting of a column of grains has been studied using a one-dimensional 

Ising model with long range directed interactions in which down and up spins rep- 
resent orientations of the grain having or not having an associated void. When the 
column is not shaken (zero "temperature") the motion becomes highly constrained 
and under most circumstances we find that the generator of the stochastic dynamics 



O ' assumes an unusual form: many eigenvalues become degenerate, but the associated 

multi-dimensional invariant spaces have but a single eigenvector. There is no spec- 
tral expansion and a Jordan form must be used. Many properties of the dynamics 

^ , are established here analytically; some are not. General issues associated with the 

^«0 ' Jordan form are also taken up. 

oo ■ 
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I. INTRODUCTION 

X- 

^ . Zero temperature dynamics is parsimonious. It is dominated by its constraints and lends 

itself to non-generic behavior. Our focus is on a model of granular compaction and we will see 
that the peculiarities of the allowed motions induce non-exponential relaxation even when 
all eigenvalues for the stochastic dynamics are real. The mathematical mechanism behind 
this uncharacteristic behavior for a Markov process lies in the fact that the associated 
eigenvectors do not span the entire state space and the best one can do is to represent 
the generator of the stochastic dynamics as a Jordan form. Although it is rare for this 
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mathematical construct to appear in a physical context, it is not unheard of [1-7], and in 
particular, in some of these there are also instances of non-exponential decay despite the 
absence of memory. 

Granular materials can exhibit features significantly different from those of traditional 
fluids [8]. A model of such materials was developed in [9-13] and also includes (strictly) 
positive temperature, which in this context relates to random shaking of the material. The 
model [9, 10] consists of a finite column of N grains, labeled by their depth n = 1, . . . ,N. 
Each grain has an orientation variable (T„ = ±1. Grain n is called up or ordered when 
0"„ = +1 and down or disordered when o"„ = —1. An ordered grain occupies one space unit. 
A disordered one traps a void and wastes space; it is said to occupy space e, so that it 
traps a void of size {1 — e) alongside it. The quantity e is effectively a shape parameter; see 
[8, 13]. Within the column we define a local field, hn, whose purpose is to model the effect of 
compacting constraints. In the simplest case [9, 10], this is only due to grains above grain-n: 

n-i f ., _ ^ 

hn = / fWm) = ^nQ — 111 , with the "shape factor" f(a) = < ' . (1) 

^ \ e, a = -1 

In Eq. (1), ni is the number of plus spins above spin-n, and no the number of minus spins 
above it. Under the dynamics (to be specified in a moment) spins tend to order in the field 
direction at finite temperature and must orient along the field at zero temperature. This 
models the observed local compaction in granular materials: the system tends to eliminate 
its voids [14, 15]. 

The system undergoes continuous Markov dynamics generated by a matrix of transition 
probabilities. In one unit of physical time the system averages A^ steps. Zero-temperature 
dynamics tends to retrieve ground states, and is defined as [9, 10]: 

an -> sign hn ■ (2) 

For vanishing hn, the simplest version of the dynamics gives a 50% chance of a switch in 
o"„. Surprisingly perhaps, all 2^ configurations can be reached, implying irreducibility of 
the transition matrix and strict positivity of the stationary state. 

Although our preliminary work using the observable representation [16-18] has proved 
useful for non-zero temperature dynamics in this model, it turns out that there is a seri- 
ous obstacle to its application when the constraints of zero temperature are imposed: the 
stochastic matrix generated under these rules does not have a full complement of eigenvec- 
tors. 

In this article we will develop in detail the stochastic matrix governing this process and 
demonstrate the aforesaid properties. Much of the power of our proofs arises from the use of 
a convenient basis. The eigenvalue spectrum turns out to be extremely simple and is shared 
by other models. We will show that appropriate combinatorial coefficients characterize the 
dimension of the space associated with each eigenvalue. Moreover, in some cases there is 
even a finer structure when, modifying the rule just given, the h = rate is changed from 
I to I -|- (5 for (strictly) positive or negative 6 (but with \6\ < |). Here too the dimension 
associated with each eigenvalue is given by a combinatorial coefficient. 

Nevertheless, there is still a great deal that is unknown in this system. For example, we 
can characterize basis vectors for only half the dimensions of the various invariant spaces 
(the construct that replaces eigenvectors). This is related to what we consider to be one 
of the important questions raised here. In quantum mechanics symmetries usually relate 



degenerate states. What is it that unites states in the same invariant space? Another 
question, not unique to this work, is, why the Jordan form? It is non-generic; in fact we 
will display results showing that even at zero temperature a slight change — in particular 
the non-zero 6 modification mentioned above, but at different e values — returns us to the 
conventional world of successful spectral decomposition. For both these questions we will 
present intuitions backed by mathematics, but — alas — not complete proofs. For example, 
the invariant spaces turn out to have dimensions given by combinatorial coefficients, and 
these will be shown to arise because of the ways of choosing a given number of spins out 
of the entire collection. So the dimension is accounted for, but there are details of the 
characterization of the space that are missing. Similarly, the Jordan form will reflect a 
cascade process in decay, as the system finds its way to the stationary state in the face of 
the many constraints imposed by zero-temperature dynamics. We are able to prove that a 
Jordan form is needed, i.e., that the eigenvectors cannot span. But as to showing that each 
invariant space has but a single eigenvector, we come tantalizingly close, but gaps remain. 

In Sec. II we present the matrix generator of the stochastic dynamics that implements the 
rules just given. With appropriate numbering of its states it satisfies a recursion relation as 
the number of spins is increased. That recursion is sufficient to deduce the spectrum of the 
dynamics and to allow generalization to other similarly structured matrices. The multiplicity 
of invariant spaces is also established. In this section we also introduce a reference model, 
a random walk on the edges of a hypercube, that plays a role later in our development. 
In Sec. Ill we examine the eigenvectors — and non-eigenvectors — associated with some of 
the eigenvalues. Following that, in Sec. IV the anomalous time- dependence is explored. In 
Sec. V we approach the dynamics from a different perspective and study the behavior of 
correlations, which because of the quasi-Boolean nature of our variables is equivalent to the 
analysis of probabilities. Using the independence of spin-A; dynamics on that of spins below 
it, we establish that for the correlations the dynamical equations assume a triangular form. 
This provides another way of discerning the spectrum and other properties, and gives insight 
into the physical basis of the properties derived from the abstract algebraic approach. In 
Sec. VI we work through the smallest non-trivial example to illustrate our general results. 

It also turns out that the zero-temperature dynamics does not always lead to a Jordan 
form. In particular, when the field h can vanish in the interior (possible for particular values 
of the parameter e) and for 6^0, the generator of the stochastic dynamics acquires complex 
spectrum, with a full set of eigenf unctions. This material is in Sec. VII. Following that, in 
Sec. VIII we take another perspective and deduce further dynamical properties based on the 
existence of a sum rule for decay rates. 

Finally, in Sec. IX we show that the basis for diagonalizing the hypercube random walk 
brings the stochastic dynamics matrix to triangular form thereby unifying the two ap- 
proaches, that of correlations and that of algebraic recursions. As for the other recursion 
results, this triangularity is more general than the particular physical model from which we 
began. 

The last section reviews and discusses our results. 



II. THE STOCHASTIC DYNAMICS 

A. Defining and labeling the generating matrix 

The central object of study is the generator of the stochastic dynamics. For convenience 
we change state-label conventions from that in Refs. [10, 11, 13]. Let the state of spin k 
counting from the top be fik = {ck + l)/2. The binary string of O's and I's corresponding 
to an A^-spin state must be ordered when the transition probabilities are written in matrix 
form. For reasons that will become evident we label states in reverse binary order. That is 
i^ = 1 + J2k=i l^k^^^^i with yu = (/ii, . . . , jim) and /xi the top spin. Here are two examples, 
iV = 2 and iV = 3: 
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(3) 



Following the rules outlined above, we give the transition probability for going from state- 
yU to state-z/. This rate will be called w{v, fi) and is the probability, per unit microscopic time, 
for the transition u ^ ^ (so one reads from right to left). It is the continuous time generator 
and its diagonal is adjusted so that "^^wlu, n) = 0. The fact that column sums of w add 
to zero already implies that it has the left eigenvector Ao(z/) = 1, with eigenvalue 0. The 
corresponding right eigenvector is the stationary state, which by virtue of the irreducibility 
mentioned earlier is non-degenerate and strictly positive [19]. By a variation of Perron- 
Frobenius theory, it is known that all eigenvalues of w have or negative real parts. Entries 
in w will be referred to either by giving n or i^ as defined above. Time evolution acts to 
the right on probability distributions ("p") in the following way: p{t) = exp{wt)p{0). 

For fi ^ u, the 2^x2^ matrix w{v^^) (or WNii^jf^) when we wish to emphasize its N- 
dependence) is zero unless the states fi and u differ in precisely one spin entry. Let that 
spin be the m^^ (from the top). Let the "change" in going from fi to u he defined as the 
binary value of the target (z/^) in the m^^ site minus the binary value of the source (/im) at 
that site. Recall that the field at level n is defined as /i„ = euq — rii, with no the number of 
down spins (strictly) above spin n and ni the number of up spins, e is a parameter. Note 
that the labels and 1 on the n's in the definition of h now correspond to the /i„ value. 
If the field and the "change" are both positive, or both negative, the transition can take 
place and WAr(z/, /x) = 1, otherwise not (in which case WAr(z/, /x) = 0). If the field is zero, 
WAr(z/, /i) = 1/2 + 6, irrespective of the change, with 6 a second parameter. 

For any particular A^, the set of e values breaks into two classes: "generic" and "non- 
generic" or "special" [10, 20]. "Generic" means the only place the field h can be zero is 
above the first spin. This corresponds to irrational e or to rational numbers that (in reduced 
form) involve sufficiently large integers. "Non-generic" e allows the field to vanish within 
the column and leads to enhanced fluctuations and to significant 5-dependent features in 
the spectrum for non-zero 6. Because in our zero-temperature model being exactly zero is 



different from being almost zero, w is stepwise constant as a function of e and its value at 
"special" e's is not its limit as e approaches these values. 



B. Recursion and eigenvalue spectrum 

The most important features of wn arise from its recursive structure as a function of A^, 
most evident in the state enumeration listed above. We can start with N = 1, 6 = 0, for 
which our assertion that 

"'^(14)' '"' 

requires a one sentence justification: the field, h, is zero, putting |'s on both off- diagonals, 
and the diagonal is adjusted to give zero column sums. For N = 2 the uppermost spin again 
sees zero field so that the upper-left and lower- right blocks are the same (i.e., as in Eq. (4)), 
except for the diagonal. This is a consequence of the numbering scheme. In particular the 
/ii values have the same pattern in both blocks; they differ only in their /i2 value. 

The off-diagonal blocks (both 2-by-2) deal with transitions in which only the deepest 
spin (in this case number 2) is changed. Hence it can only have entries on its diagonal. 
These entries are 0, 1 or 1/2 depending on the sign of the field or whether it's zero. Before 
correcting the diagonal for zero column sum, W2 has the following appearance: 

. . I \ *^1 y non-diagonal portion ^— ^ 1 / 1^ \ 

\ *^2y non-diagonal portion I 'T' , ■. I 7 \ / 

i_i \ *^1 y non-diagonal portion / 

where A and A are themselves diagonal matrices (of the same size as wi). What is important 
to note is that A + A = 1, the size-wi identity matrix. This is because for each pair of states 
on A's diagonal (which differ in a single spin) if the field is positive for A it is negative for 
A so one matrix element is 1, the other zero. If the field is zero, it is zero for both, and both 
are 1/2. A moment's reflection shows that getting the column sums right is also easy and 
the form of W2 is 

""'=[ A w,-a)^ (') 

The exact form of A — which is to say, where it has zeros, ones and halves — will depend on e. 
(In Sec. VII we take up the non-zero 6 case, which for non-generic e can violate A + A = 1.) 
The arguments we have just given for building W2 from wi are valid for any A^ and by 
the same reasoning one has 

I wn - An An . /„v 

y An wn- An ' 

Here we have added the label A^ to A indicating that it is a 2^-by-2^ matrix, as is wn- (A 
word of caution on notation: An enters the (A^ + l)-spin dynamics; it is not part of Wn-) 
Moreover, 

An + An = In, (8) 



for S = OT for non-zero S in the case of generic e, since the field cannot vanish within the 
column. (Note that 1^ is not the N-hj-N identity, but the 2^-by-2'^ identity. Similarly Otv 
is the 2^-by-2^ matrix of zeros.) 

It is remarkable that the full eigenvalue spectrum as well as the dimensions of the invariant 
spaces (whether a Jordan form is needed or not) can be deduced from this recursion alone 
(plus properties of Wi that allow an induction). The argument proceeds by examining the 
characteristic polynomial of wn- Let 



Pn{\) = det (wn - Al 
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The step from Eq. (12) to Eq. (13) depends on Eq. (8). (As indicated, this holds for generic 
e and any 6, or for non-generic e and zero S. It does not hold for non-generic e and non-zero 
S. For that case see Sec. VII.) 
The inductive hypothesis is 
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(15) 



fc=0 



where C^ = N\/k\{N — k)\, the combinatorial coefficient. This assertion is trivial for A^ = 1 
or 2. 

From Eq. (14) it follows that 
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The inductive hypothesis is now proved by observing that 
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It follows that the eigenvalues are {0, —1, . . . , — A^}, with eigenvalue —k having an invariant 
space of multiplicity C^ . 



Because of our analytic information on the spectrum and multiplicities, this situation will 
be referred to as "integrable." This is to be contrasted with the non-generic e, non-zero S 
case, which we will characterize as "non- integrable." 

Remark 1: The spectrum and invariant space structure are the same for any family of 

matrices obeying Eqs. (7) and (8) for which the induction can be initiated. An illustrative 
example is Aj. = ^Ifc for k < N — 1. Call the associated matrix y^. This is the generator 
of a random walk on the A^-cube. It can also be looked upon as the infinite temperature 
limit of an A^-spin Ising model in which only single spin flips are allowed. Unlike our zero- 
temperature granular dynamics, the matrix y is symmetric and has a full complement of 
eigenvectors. They are given by the following construction. Let A be a subset of the numbers 
1 through A^. There are 2^ such subsets. For this subset define a particular 2^-vector, va, 
as follows. Let jj, = {ni, . . . , (In) be an A^-string of O's and I's. Then fyi(/i) = (— 1)"(/^'^), 
where n(/i. A) is the number of I's in /i that fall in the subset A. (In other words, if Xa i^ 
the characteristic function of A, n{fi,A) = J2kf^k ■ XaW-) ^^ i^ '^ot difficult to show that 
this is indeed an eigenvector of yN and has eigenvalue —\A\, with vertical bars indicating 
cardinality. This immediately implies that the multiplicity of the eigenvalue —k is C^ . A 
propos the A^-cube interpretation, when the image of this random walk under the observable 
representation [21] is plotted, one gets in fact a cube, not at all surprising in view of [18]. 

Remark 2: Many other matrices having the same eigenvalues and invariant space dimen- 
sions can be constructed by varying A^ as you build toward A^. Some require the Jordan 
form, some do not. There is no need to have the entries in A be confined to {0, 1/2, 1}, nor 
even to the reals. There is no need to start from the wi given above. For example one can 
start from 



2. Wi = [ ). This has eigenvalues and 1. It is only stochastic for 



1. Wo = 0, a 1-by-l matrix. It has the single eigenvalue 0. 

a fj,{a + 1) 

-a/fi —(a + 1) 
/i = 1. 

So there is also no requirement that the matrix be stochastic. The matrix A can be anything 
of appropriate size. Possibility 2 is more general than 1. The induction can also be started 
at larger matrices, say, W2, allowing yet larger classes of isospectral operators. Approach 2 
reduces to Approach 1 for /i = l,a = — 1/2. 

Remark 3: In our numerical experience, all N > 1 matrices constructed under the rules 
for granular dynamics were maximally Jordan, by which we mean that there are degenerate 
eigenvalues and that each invariant space has but a single eigenvector. (This includes generic 
e, (5 7^ 0, but does not apply to the non-integrable situation, i.e., non-generic e, non-zero S, 
where the eigenvalues cease to be degenerate.) 

Remark 4: The same operations by which we step-by-step simplified the determinant 

can be performed on the original matrix, casting additional light on its structure. See 
Appendix A. 

Remark 5: In building matrices using the recursion Eq. (7), if the entries of the diagonal 
matrix A are taken randomly from the set {0, |, 1}, one often has non-spanning eigenvectors, 
i.e., the need for Jordan forms. However, the invariant spaces are not necessarily maximal, 
i.e., there can be more than one eigenvector with given eigenvalue. We further remark that 
discerning whether or not a matrix requires a Jordan form can be numerically delicate. One 



indicator is that when one attempts to diagonahze by conventional methods one finds that 
the associated "eigenvalues" contradict analytically determined properties, for example, in 
having a non-zero imaginary part (see Sec. VII A for the explanation). 

We next consider non-zero 6, generic e, so that the only site that can experience zero-field 
is the top one. The only transitions affected are those between an odd-numbered site and 
the even-numbered site with index one larger. Then the change in w due to the presence of 
6 is 

/a ...\ 
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\0 ... aj 



with a = ( ' J . (19) 



There are 2^ ^ copies of a in wn- Since w is linear in 5, WNiS) = WAr(O) -|- S^^. For 
Pn{\,S) (defined as det {wn{S) — \In)) we make the inductive hypothesis: 
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-,N- 



{X + kf^ {X + k + 26f'^-' . (20) 

fc=0 

(Recall (Eq. (18)) that Cj^ = C^^^ + C^Zi , so for 5 = this reduces to the former case.) 
This hypothesis can be directly verified for A^ = 2 or 3, and follows from our Eqs. (10)-(14) 
(which, as indicated, remain true for generic e) together with identities of the form Eq. (18). 
In Sec. VG we provide another way of reaching the same conclusions. 

It follows that the invariant subspaces are now smaller: new eigenvectors emerge, one per 
invariant space. The span of the (5 = invariant space associated with the eigenvalue —k is 
not the sum of those for eigenvalues —k and —k — 25 for non-zero 5. 

III. EIGENVECTORS AND INVARIANT SPACES 

To address the nature of the eigenvectors and invariant spaces additional tools will be 
used. It will also be useful to work with the transpose of w, since its right eigenvectors often 
have simpler structure. Let qn = wj; (with possible suppression of the index A^). 

We first observe that there is a doubling rule as A^ increases. Suppose u^ is the (true) 
eigenvector of qn of eigenvalue A. Then it is easy to show using Eq. (7) that 

In other words, the doubled true eigenvector is a true eigenvector of the next larger (by a 
factor 2) "^." 

Remark 6: For w the doubling is slightly different. There is a minus sign and a shift in 
eigenvalue. One can immediately verify that if w^v^ = Xv^ , then 

^^+1 ( 5 ) = (^ - 1) f 5V (22) 



But the doubling is more general and applies to the entire invariant space. Let -uj^^ be an 
element of the invariant space of A (for the matrix g), with the convention that i = 1 is the 
true eigenvector, so that 1 < i < Nx, with Nx the dimension of the invariant space (which 
is Ci^i for 6 = 0). The others are labeled according to their place in the Jordan form in the 
following way: 

((7JV - A/^X, = <,„! , for2<£<iVA. (23) 

Note that in writing u^^ we will often replace A by its negative, since this is unambiguous 
and allows the integer structure to stand out. Now consider Qn+i applied to the doubled 
vector (for i >2), 
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(24) 



The final result is still a doubled vector, as a result of which the A and A entries in g^+i 
have no effect on it. Therefore one can continue to apply {gN+i — ^^n+i), each time reducing 
i, until, reaching 1, the vector is annihilated. Since this vector was in the invariant space 
for N it will be annihilated in at most C^i steps and therefore is certainly in the invariant 
space for A^ + 1. Physically, the feature arises because the addition of another spin does 
not modify the transitions between spins that are above it [10]. As one goes from N to 
N + 1 spins, the spin states are doubled; the invariant space for a particular eigenvalue of 
the iV-spin configuration is also doubled to reflect this in the A^ + 1-spin transition matrix, 
but as all the corresponding transitions are unaffected by the presence of spin A^ + 1, it is 
left otherwise untouched. 

On a lighter note, the doubling property means that for each A^ + 1, 2^ dimensions (out 
of 2^"*"^) are accounted for, so you might say that half the job is done. Unfortunately there 
are many interesting questions in the other half. 

Remark 7: The operations and relations discussed here can be conveniently phrased in 
terms of the operators defined in App. B. 

As observed earlier (Remarks 1, 2 and 5), the recursion, Eq. (7), establishes the spectrum 
and multiplicity of the invariant space, but it does not indicate the number of true eigen- 
vectors in each invariant space. For example, in the symmetric case (A = lAr/2) there is a 
full complement of eigenvectors. We have found however, that all invariant spaces arising 
from the zero-temperature granular dynamics model have but a single eigenvector (also for 
(5 7^ 0, e generic). We know this numerically up to A^ = 6. We next show analytically that 
for our model the invariant spaces do contain non-eigenvectors (a statement that is much 
weaker than our numerical experience). 

For A^ = 2, A = —1, one obtains from the Jordan form representation that (uli)^ = 
[1,-1,1,-1] and (^1^2)^ = [2,0,0,-2]. It follows that uf-^ = (wj + l) -uj 2 ^^^ that UI2 



10 

is not an eigenvector. Now we look at the A = — 1 subspace for A^ = 3 and study the 
doubhng of m^ 2- ^^is corresponds to the separate action of wj on each portion. This does 
not produce zero, but instead produces the A^ = 3 eigenvector for eigenvalue —1; i.e., the 
doubled u1 2 is in the invariant space, but is not an eigenvector. This argument obviously 
continues to hold for higher A^. Note that the doubled objects we have produced may not 
coincide with what appears when obtaining the Jordan form. This is because arbitrary 
pieces of (e.g.) the eigenvector can be added without changing the lowering property (as 
exemplified in Eq. (23)). 

A. Largest magnitude eigenvalues 

The simplest eigenvector of g is that with eigenvalue 0, namely the vector whose entries 
are all 1. There is corresponding explicitness for the largest magnitude eigenvector, but 
this time we look at the eigenvector of w, not its transpose. Recall that for wn the largest 
magnitude eigenvalue is A = — A^. 

Suppose then that 

=.««(;)=a(;). (25) 

with a and f3 column vectors of length 2^. Then from the recursion, Eq. (7), it follows that 

w^a - Aa + Af3 = \a (26) 

wn/3 - A/3 + Aa = A/3. (27) 

Adding these one gets 

WN{a + P) = X{a + P) . (28) 

Now suppose A = — A^ — 1. We know that — (A^ + 1) is not in the spectrum of w^, implying 
that P = —a. Now for wi of Approach 2 (of Sec. JIB) the eigenvector with eigenvalue —1 is 

irrespective of the value of a in that matrix (our case is a = —1/2). It follows by induction 
that for the largest absolute value eigenvalue of wn the eigenvector consists entirely of ±/i's 
and ±l's. The pattern is built by successive attachments. For /i = 1 (which is the value 
for zero-temperature granular dynamics) the entry is (— 1)"^^ with z/i, the number of I's in 
the state. For /i 7^ 1, the odd entries are multiplied by fi. If one uses spin notation, i.e., 
fi£ ^ (X£ = 2fi£ — 1, then this largest magnitude eigenvector is given by the product of the 

spins, 0'i(T2 • • • 0"Ar. 

The simplicity of this eigenvector is to be compared to the corresponding eigenvector of 
gN- First note that by appropriate multiplication by powers of 2, all vectors in the invariant 
spaces can be written as integers. We have observed that (with minimal multiplication) 
the elements of these vectors grow rapidly. For example, the first component of Uq ^ is 
2075566815213212256361. Factoring this number (it equals 19 x 41 x 193 x 239 x 18523 x 
3118406479) or others like it has not given us clues to its structure. The algorithm for 
obtaining eigenvectors involves sums and products of the matrix elements, so it would appear 
that the many halves in gi\i are not canceling in any systematic way. It should be noted that 
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this growth apphes not only to ■uj^^(l), but to pretty much any component for which there 
was no neat exphcit form. 

Remark 8: Consider, for wn+i, the eigenvalue N, which is an eigenvalue of wn as well. 
Then Eq. (28) implies that a + /3 is an eigenvector of wn- Assuming that a + /3 is not zero, 
then /3 = u^ i — «, with u^ ^ a true eigenvector of W]\[. 

B. On the "maximal Jordan" property 

All zero-temperature granular models give rise either to what we have called "maximally 
Jordan" spectrum or to complex spectrum. Thus all e and 6 except those combinations 
yielding fish bone spectrum (i.e,. generic e or non-generic e, zero S), are maximally Jordan. 
(Recall, this means that degenerate eigenvalues, of which we assume at least one present, 
have but a single eigenvector. Recall too that \6\ < 1/2.) This assertion is based on numerical 
solutions up to and including A^ = 6. Included in this observation are the smaller subspaces 
for e generic, 6 nonzero. 

On the other hand, if we use the recursion of Eq. (7) and generate A's whose diagonal 
entries are O's 1/2 's and I's (but which need not be realizations of the granular model), 
sometimes one gets the maximal Jordan property, sometimes not. In this subsection we 
explore further criteria for this property. 

Consider the eigenvalue equation for wn+i (A^ + 1 spins) 

.wn-An An \ f x\ _ f WNX-X + AN{x + y)\ _ , f x\ 

with X and y 2^-vectors, k = 0,1, ... ,N and use has been made of A + A = 1. We only 
consider 6 = 0. Add these equations to get 

WN{x + y) = -k{x + y) . (31) 

We make the inductive hypothesis that every invariant space of w has hut a single eigenvec- 
tor, which in the presence of degenerate eigenvalues (which we know occur for A^ > 1) is the 
maximal Jordan condition. Then Eq. (31) implies that x + y = aVk, with v^ the eigenvector 
of Wn of eigenvalue —k. a is a real number and there are two possibilities: it is zero or it is 
not. If it's not, renormalize x and y to make it —1. 
Case I. a = 0. In this case x = —y. From Eq. (30) we have 

wnx = —{k — l)x , (32) 

so that X is an eigenvector with eigenvalue —k + 1 and we have (re-) discovered the doubling 
property, namely that once you have a true eigenvector of eigenvalue —k you can get the 
(or "a" if it's not maximally Jordan) eigenvector (of eigenvalue —{k — 1)) for one more spin 
by doubling and changing sign. 
Case II. a = —1. Then x + y = —v^ and we use Eq. (31) to obtain 

wnx + {k — l)x = ANVk 

WNy + {k - l)y = ANVk (33) 



12 

To prove the maximal Jordan property we need to show that Eq. (33) has no solution. In 
that way, moving to higher N, you never get a new eigenvector, just doubling of the old 
ones — with one exception. That exception is k = and that's the new eigenvector (of w, 
not w^). 

Remark 9: For k = 0, Eq. (33) can always be solved, since (+1) is not an eigenvalue 

of w, and the left hand side is invertible. As remarked elsewhere in this article, the k = 
eigenvector of w, when scaled to involve integers with no common denominator, has entries 
that grow faster than exponentially (and to us, unpredictably) with N. These entries, 
normalized to sum to unity, are the probabilities of the various spin configurations in the 
stationary state. As indicated below, the nature of the stationary state has a sensitive 
dependence on the parameter e and on the boundary conditions. 

Considering now only the k > case, the left hand side of Eq. (33) has an important 
property: it cannot produce a vector proportional to ffc-i,n, where n = C|^i, and we use the 
notation (of Sec. Ill) indicating that f fc-i,„ is the deepest vector in the Jordan chain. This 
vector is characterized by the property that it, and only it, can survive n — 1 applications 
of {w + {k — 1)1). (It is not unique, but can be selected uniquely by the demand that it 
be orthogonal to the rest of its invariant space. But even without selection the foregoing 
criterion obtains.) 

So the problem of showing this system to be maximally Jordan reduces to showing that 
either Avk or Avk has non zero overlap with Vk-i^n (with n = C^_i), the deepest vector in 
the Jordan chain.. 

Remark 10: This condition is necessary and sufficient. This means that all the maximally 
Jordan matrices generated with A's that do not correspond to what you'd get from the 
granular model also have no solution to the above equation. So the criterion is more general 
than the model, but for dynamical reasons the model satisfies it with each additional spin. 

Remark 11: Like the stationary state, this lowest vector in the Jordan chain is created 
de novo with each succeeding increase in A^, and is not obtained by doubling. This is clear 
from Eq. (24), since vectors in the invariant space obtained by doubling cannot require the 
maximal number of steps for the increased A^ to reach the eigenstate under applications of 

(w + kl). 

Remark 12: A word of caution regarding the establishing of whether or not Ax or Ay 
have components along some direction in the vector space: as usual when the Jordan form 
is needed, many of the usual tools break down. Thus even for a general stochastic matrix, 
if it does not require a Jordan form you can still have leftj, ■ right^ = 5ai3, where "left" 
and "right" are eigenfunctions of w and w^ transpose respectively (so that "left^" is a left 
eigenvector of w). When w requires a Jordan form, this is not true and for our maximally 
Jordan matrices the left and right eigenvectors (when the invariant space is of dimension 
greater than one) are orthogonal. 

Remark 13: Because the doubling properties of w and its transpose ("(?") differ slightly, 
the corresponding equations for g take the following form 

gNX + kx = A^Uk-i 

gwy + ky = -A^Uk-i. (34) 

The condition for proving the maximal Jordan property is that application of A or A to the 
next slower eigenvector (i.e., —A = A; — 1) not yield any component along the bottom of the 
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Jordan chain for eigenvalue —k. This condition is particularly simple for k = 1 (since Uq is 
all I's), and we have noted numerically that the lowest in the k = 1 Jordan chain is indeed 
highly correlated with the associated A. 

IV. TIME DEPENDENCE 



Non-monotonic time evolution occurs in this model, just as for a critically damped oscil- 
lator. From the definition of w, the time dependence of an A^-spin system is 



p{t) = exp{wNt)p{0) 



(35) 



with p the (vector) probability distribution (at times and t). Were this a situation with 
(only) bona fide eigenvectors one would expect p{t) to be a sum of terms of the form 
exp(Aat) with a running over the eigenvalues. With eigenvalue degeneracy, terms of the 
form t'' exp{\at) can occur, for k up to one less than the level of degeneracy. In this section 
we show precisely how the power law enters in the probability distribution. In Sec. VI it 
will emerge naturally in a study of correlations. 

Let the initial state be of the form po + pvf^ with po the stationary distribution (a.k.a. 
f(^i) and p small enough so that all components of p(0) are non-negative. The action of w on 
this is as follows: it annihilates poJ the vf^ term is also annihilated but two new terms replace 
it: one is proportional to the eigenvector (v^-^) and the other is vf^_^. The exponentiation 
of w is not quite so simple as when there is a spectral decomposition, but it remains true 
that this process respects the invariant spaces. Moreover, one can immediately verify that 
for 



m 



/A 1 

A 1 

... 

VO ... 






A 




\ 



1 



(36) 



the matrix exponential is 



exp(tm) = e 



t\ 



( 1 t e/2 tV6 

1 t tV2 

... 1 

yo ... 



t^-V(£-l)!\ 



(37) 



/ 



*'^ exp(tm) is upper triangular with the k^^ ascending diagonal constant and equal 
= 0, ...,£ — 1, with I the dimension of m. To see the power law for the initial state 



That is, e 

totV^!,^ 

given above, we require some notation. Let the Jordan decomposition of w be accomplished 

by the following transformation 

wV = Vj , (38) 

with j composed of blocks of the form given in Eq. (36). Let u^ be the (row) vector satisfying 
1 and which is orthogonal to all other vectors (it is given by an appropriately 



T„.N 



U V 



e,i 



normalized row of V ^). With this notation, the object having power law decay is s{t) 
u^ exp{wt)v^i. This is illustrated in Fig 1. 
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FIG. 1: Survival amplitude, s{t), of an initial state of the form, po + pvf 2- The stationary state po is 
removed and an appropriate projection of the remainder is taken. What is shown is the logarithm 
of what's left, both in raw form and after division by t^. Note that this is the 3-spin generating 
matrix and we are looking at the invariant space associated with the eigenvalue -1. 

V. DYNAMICAL EQUATIONS FOR SPIN CORRELATIONS 

A. Background 



In this section we investigate the dynamics of the model based on (equal-time) spin 
correlations. These quantities are defined as the time-dependent mean values of products of 
spins: 

M,,..,,(t) = K(t)---a,,(t)), (39) 

where ii, . . . ,ik is an ordered A;-tuple of distinct labels in 1, . . . , iV. It is sufficient to consider 
multilinear functions, i.e., functions which are at most linear in each of the spins, since each 
spin obeys af = 1. There are C^ such correlations at level k, and hence 2^ correlations 
in total (including M^ = 1). The spin correlations thus form a basis of observables. The 
knowledge of those correlations is tantamount to that of all configuration probabilities. 

The temporal evolution of spin correlations is dictated by first-order linear differential 
equations, which can be written down explicitly. 

To illustrate the approach, we first consider the infinite-temperature situation, i.e., the 
random walk on the hypercube in dimension N, mentioned in Remark 1. In this case, 
each spin fiips at Poissonian times with unit rate, independent of the others. Consider for 
definiteness the evolution of the first spin cri(t) during an infinitesimal time interval dt. We 
have 

/, i,\ f <^iii) with prob. 1 — dt, ,,^. 

^'(' + *' = |-.,W with prob, it («' 
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The average, Mi(t) = (ai), therefore obeys Mi(t + dt) = (1 — 2dt)Mi(t). As a consequence 

dMi 

"dT 

and more generally 

dM 



-2Mi, (41) 



^P = -2A;M,...,,. (42) 

The basis of spin correlations therefore diagonalizes the dynamics. We thus readily obtain 
the spectrum of the Markov matrix. The eigenvalue A = —2k has multiplicity C^ (in 
accordance with the number of possible subsets associated with Eq. (42)), and the quantities 
^n..ife provide an explicit set of eigenvectors for the infinite temperature model. This 
diagonalization is equivalent to that described in Remark 1. 



B. General approach 

Let us apply the above approach to our column model. 

• If /i„ 7^ 0, (T„ is updated with unit rate according to the rule 

(T„ -> sign/i„. (43) 

• If hn = 0, (T„ is updated with rate w = j + 6 [22] according to the rule 

(Tn -^ -(Tn- (44) 

For definiteness, we start with the following three hypotheses: 
(1) Free boundary conditions (i.e., /ii = 0). 
{2) w = \ (i.e., (5 = 0). 
(3) e is generic (i.e., irrational or 'large' rational). 

As discussed earlier, the phrase "e is not generic" means that e is one of the special 
values that allows vanishing h below the first spin (see Table I), i.e., for given system size, 
N, e = p/q (irreducible) with p + q < N — 1. 

For w = |, the rule (43) still holds on average for /i„ = (defining signO = 0). Flipping 
spin an with probability ^ sets a^ = ±1 with equal probability. Therefore in full generality 

_^ = _M„ + (sign/i„) , (45) 

where, as above, M„ = {cXn). The same observation will allow us to write down differential 
equations for all the observables, as they are multilinear functions, i.e., linear functions of 
each of the spins. 

The above approach is to be put in perspective with the pioneering work of Glauber on the 
dynamics of the ferromagnetic Ising chain [23]. Although Glauber only considered one-spin 
correlations (magnetization profile) and two-spin correlations in his original work, several 
subsequent papers have been devoted to a systematic extension to an arbitrary number of 
spins [24-27]. 

Before considering the general situation, it is useful to first derive explicit equations for 
the first few values of N. 



16 



C. The case iV = 1 

There is only one non-trivial correlation, Mi(t) = (o'i(t)), which obeys the differential 
equation 

_-i = -Mi + (sign/ii), (46) 

with 

sign hi = 0, (47) 

so that 

This equation for Mi still holds for larger system sizes, because the motion of spin-1 does 
not depend on those below it. We thus obtain 

Mi(t) = Mi(0)e-*. (49) 



D. The case iV = 2 

Besides Mi{t), there are two new correlation functions at size N = 2, M2{t) = {o'2{t)) 
and Mi2{t) = {ai{t)(72{t)). They obey the differential equations 

dMa ,, , . , , 

— -— = -M2 + (sign /i2) , 
dt 

—^ = -2Mi2 + (sign hi-a2) + {ai ■ sign /is) , (50) 

with 

sign/i2 = -o-i . (51) 

Note that this holds irrespective of e. It follows that 

dMs 



dt 
dMi2 

dt 



-M2 - Ml, 

-2Mi2 - 1. (52) 



These equations continue to hold for larger system sizes, because of the dynamics of a given 
spin is independent of those below it. 

A detailed analysis of the dynamics of the model with N = 2, albeit with a general 
w = ^ + 6, will be performed in Section VI. 



17 



E. The case iV = 3 

There are 4 new correlation functions at size A^ = 3, which obey the differential equation 
dMg 



dt 
dMi3 



-M3 + (sign h) , 

-2Mi3 + (sign hi ■ 0-3) + {ai ■ sign h^) , 



dt 

-2M23 + (sign /i2 ■ 0-3) + (0-2 • sign h^^) , 



dM23 



dt 

dMi23 

dt 



with 



-3Mi23 + (sign hi ■ (J2(Tz) + {(Ji ■ sign /i2 ■ 0-3) + {aia2 ■ sign /i3) , (53) 



sign/13 = -(0-1 0-2 - 0-1 - 0-2 - 1). (54) 



This identity holds whenever e <1. Assuming this for definiteness, we have 

^ = -M3 + i(Mi2 -M1-M2- 1), 
dt 2 

' 2Mi3 + -(-M12 - Ml + M2 - 1), 



dt 2 

dM23 .,, ,, 1 



-2M23 - Mi3 + -(-M12 + Ml - M2 - 1), 
dt 2 

^^ = -3Mi23 - M3 + -(-M12 - Ml - M2 + 1). (55) 

F. The general situation 

The structure of the dynamical equations is apparent from the above examples. The 
general procedure consists of the following steps. 

A general spin correlation Mj^..,,^ obeys a differential equation of the form 

where 

H\l\.i^ = (ffii • • • f^i,_i ■ sign hi. ■ ai^^, .--(Ji^) (57) 

is obtained from Mi^,,,i^ by replacing the j^^ spin a^. by sign/ij,. 

For each depth n, sign/i„ is a symmetric multilinear function of the spins above n, i.e., 
o"i, . . . ,0"„-i. Its explicit expression can be obtained by listing the values of sign/i„ for all 
the 2""^ spin configurations, and fitting this set of values to a multilinear function with 
arbitrary coefficients. The resulting expression, generalizing (51) and (54), depends on e. 
More precisely, it depends on the relative position of e with respect to the special values 
that enter at depth n, such that /i„ can vanish. These are the {n — 2) rationals e = p/q 
(not necessary irreducible) such that p + q = n — 1. The expression of sign/i„ is constant 
for e between any two consecutive special values. Table I gives an ordered list of the special 
values of e corresponding to the first few values of n. 
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n 


special values of e 


3 


1 


4 


i 2 

2' ^ 


5 


i 1 S 
3' -^' "^ 


6 


12 3/1 
4' 3' 2' ^ 


7 


i 1 1 2 'S 

5) 2' ■"■' ^' "^ 


8 


1 2 3 4 5 f? 

6' 5' 4' 3' 2' " 


9 


1 1 3 1 5 o 7 

7' 3' 5' ' 3' ' 



TABLE I: Special values of e when the total number of spins is n. 



Now, for a fixed generic e, let us replace in (57) sign hi- by its explicit expression in terms 
of the spins, carry this into (56), and expand the sum. The expression thus obtained is a 
linear combination of operators whose labels are strictly smaller than that of Mj^ j^. The 
label i of an operator (spin product) o"j^ . . . o"j^ is defined as 



)ji-i 



+ 2^2-1 + ■ ■ ■ + 2 



«fc- 



(58) 



This labeling provides a universal ordering of operators for all system sizes, whose beginning 
is given in Table II. The ordering of £ = 11 and £ = 12 is the first surprising or unexpected 
one, in the sense that the number of spins involved decreases at fixed A^. 



N 


i 


operator 








1 


1 


1 


0-1 


2 


2 


(^2 




3 


(Ji(T2 


3 


4 


0-3 




5 


CTlCTs 




6 


0-2 0-3 




7 


Cri(T20'3 


4 


8 


£74 




9 


0"lO"4 




10 


0-2 (T4 




11 


CJl(T20"4 




12 


0-3 (T4 




13 


aia^a^ 




14 


(T2Cr3(T4 




15 


(Ti(72(73(T4 



iV 


i 


operator 


5 


16 


0-5 




17 


aia^ 




18 


0"20-5 




19 


0'lO'20"5 




20 


^^30-5 




21 


0"! 0-3(75 




22 


0"20-3(T5 




23 


0"l(T2Cr3(T5 




24 


0-4175 




25 


0"lO-4Cr5 




26 


0-2 0-4 0-5 




27 


O"! 0-2 0-4 0-5 




28 


0-30"40-5 




29 


cri(T3Cr4(T5 




30 


0-20"30-4(T5 




31 


0"lO-20-30-40-5 



TABLE II: Universal labeling of operators. 



The dynamical equations (56) therefore have a triangular form, for any finite system 
size A^. The explicit form of those equations depends on e, albeit their triangular structure 
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is robust. This feature provides an alternative proof of the pattern of degeneracies of the 
Markov matrix studied earher by algebraic techniques. We indeed obtain at once from 
the triangular form of (56) that the spectrum consists of the negative integers, i.e., A = —k 
(/c = 0, . . . , A^), with combinatorial multiplicities C^ . Furthermore, although the occurrence 
of a single Jordan block of maximal length C^ at each level k is not proved by this approach, 
this generic feature is suggested by the large number of off-diagonal terms generated by 
expanding the sum in the right-hand side of (56). 



G. Extensions 

We extend the above approach by investigating whether and how the three hypotheses 
made below (44) can be lifted. 

• Lifting hypothesis (1) is easy. Let us remark that hypothesis (1) corresponds to hav- 
ing free boundary conditions. Indeed the field hi vanishes identically, so that the 
uppermost spin ai is free. A natural alternative consists of fixing the uppermost 
spin, along the lines of our previous works [9-12]. Let us rename this spin o"o and 
set (To = +1 for definiteness. This amounts to adding a constant to all the fields 
{hn -^ hn + /(o"o) = hn + /(+1) = Ki — !)• Accordingly, the system size (number of 
degrees of freedom) is reduced from A^ to A^— 1. The explicit multilinear expressions of 
sign hn in terms of the spins are therefore modified. The overall construction however 
still holds. The spectrum is left unchanged (up to the size reduction from A^ to A^ — 1). 

• Lifting hypothesis (2) {w = |) is also easy, provided hypothesis (3) is maintained. 
Indeed, as long as e is irrational, the only point of zero field where /i„ can vanish is 
n = 1. As a consequence, the zero-field rate w = j + 6 only affects the dynamics of 
the uppermost spin o"i. Therefore, it only enters the diagonal term of the differential 
equation (56) for the operators involving o"i (i.e., such that ii = 1) as follows: 

k 

H = l: ^^ = -{k + 25)M,..., + 5:<.,, 

i=i 

k 

^1^1: ^^ = -kM....^. + E Hit. ■ (59) 

The first line yields A = — (A; + 25) with multiplicity C^Zi for A; = 1, . . . , A^, whereas 
the second one yields A = —k with multiplicity C^^^ for A; = 0, . . . , A^ — 1. We thus 
recover the known split spectrum, derived algebraically from Eq. (20). 

• Lifting hypothesis (3) (e generic) is more subtle. The situation where e is rational 
and w 7^ I (i.e., 5 7^ 0) is indeed exceptional, in the sense that the spectrum is not 
integrable. The spectrum will be shown later to generally have a complex 'fish bone' 
structure. The gist of this complexity is as follows. For e = p/q (irreducible), there 
are internal points of zero field, that is, hn can vanish. This occurs at depths n such 
that n — 1 is a multiple of the period {p + q). As a consequence, for e rational and 
w 7^ |, the total exit rate from a configuration depends on the configuration. Hence 
there is no Glauber rule to construct differential equations. The whole construction 
breaks down. 
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To summarize, the approach based on spin correlations works everywhere except in the 
non-integrable situation where e is rational and w 7^ |, where there are good reasons why 
it cannot work. 

Table III recapitulates our results on the spectrum of the Markov matrix, as a function 
of e and 5 = w — |. The last column gives the exact expression of the average absolute 
eigenvalue A (see (72)), calculated later as a sum rule. 



Free boundary conditions 



e 


6 


spectrum 


A 


generic 


6 = 


X = -k mult. Cf (A; = 0, . . 


.,iV) 


N 
2 


generic 


(5/0 


1 ^ = "^ mult. Cf-i {k-- 
\x = -{k + l + 26) "^ ^ 


= 0,...,iV-l) 


f + ^ 


special 


6 = 


\ = -k mult. Cf (fc = 0, . . 


.,iV) 


2 


special 


(5/0 


complex 'fish bone' 


f +^7V(5 


Fixed boundary conditions 


e 


5 


spectrum 


A 


generic 


5 = 


\ = -k mult. Cf (fc = 0, . . 


.,iV) 


Af 
T 


generic 


5^0 


\ = -k mult. Cf (A; = 0, . . 


.,iV) 


N 
2 


special 


6 = 


\ = -k mult. Cf (A: = 0, . . 


.,iV) 




special 


6^0 


complex 'fish bone' 


f + 5iv(5 



TABLE III: Spectrum and average absolute eigenvalue A (see (72)) of the Markov matrix of a 



system of size N , as a function of e and 6 
boundary conditions. 



w 



\. Top: free boundary conditions. Bottom: fixed 



VI. FULL DYNAMICS IN THE CASE OF TWO SPINS 

In this section we give the full solution to the dynamics of the model in the first non- 
trivial case of two spins (A^ = 2). In this situation, the correspondence between the spin 
correlations and the configuration probabilities ^0-10-2 reads 



I.e., 



Me 


= P++ + P+~ + P-+ + P- 


Ml 


= P++ + p+_ - p„+ - p- 


M2 


= P++ - P+- + P-+ - P- 


M12 


= P++ - P+- - P-+ + P- 


++ = 


i(l + Mi + Mi + Mi2), 


+- = 


i(l + Mi-Mi-Mi2), 


-+ = 


i(l-Mi + Mi-Mi2), 


= 


i(l_Mi-Mi + Mi2). 



(60) 



(61) 
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At variance with the above general analysis, we consider an arbitrary zero-field rate 
w = ^ + S. Therefore, irrespective of e, the upper spin ai flips freely (see (47)) with rate w, 
whereas the lower spin o"2 flips with unit rate under the action of sign/i2 = — o"i (see (51)). 
The dynamical equations for the spin correlations therefore read 



dMi 
dM2 

~dr 

dMi2 
dt 



-(l + 2(5)Mi, 

-M2-M1, 

-2(1 + 6)Mu - 1. (62) 



These dynamical equations have the expected triangular form, and yield the expected spec- 
trum: A = 0, A = -1, A = -(1 + 26), A = -2(1 + 6). 

For arbitrary initial conditions, the general solution of the above equations reads 

Mi(t) = Mi(0)e-(^+2^)*, 
M2(t) = M2(0)e-*-Mi(0)i?(t), 
Mu{t) = Cu + (Mi2(0) - C12) e-2(i+^)*, (63) 



C12 = (^i^2),tat = -wrr-^ = -TTT-^ (64) 



where we have introduced the stationary- state correlation 

1 

2(1^ 5) 2w + l 

and the off-diagonal response function 

R{t) = ^^e-\ (65) 

The stationary-state correlation C12 is always negative, and it depends on the rate w in a 
smooth way. The regimes of small and large w can be understood as follows. In the w — )■ 
limit, o"i moves very slowly, and so (72 follows sign/i2 = — o"! almost perfectly adiabatically, 
so that o"2 ~ — o"! and C12 tends to —1. In the opposite regime {w — )■ 00), o"i moves so fast 
that o"2 hardly feels a driving fleld, so that C12 approaches 0. 

The response function R{t) is more interesting. For 5 = it reads 

R{t) = te-\ (66) 

The resonance phenomenon recognized through the presence of the 'secular' prefactor t 
takes place precisely at the point where the eigenvalues A = — 1 and A = — (1 -|- 25) become 
degenerate, so that a Jordan block of size Cf = 2 is needed, in agreement with the general 
theory. 

The response function increases from zero, reaches a maximum, and falls off to zero, 
irrespective of 5. The integrated response, 

decreases smoothly as a function of 5, just as the response function at any flxed time t. The 
resonance phenomenon at 5 = does not yield any particular attribute for the response 
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function. This point can however be shown to demarcate between two different regimes of 
asymptotic decay: R(t) indeed falls off as e~* all over the range S > 0, i.e., w > ^, whereas 



its decay is slower, as e^ 



-(l+2(5)t 



-2wt 



, for -i < (5 < 0, i. 



2' 

e., < w < |. These features are 



illustrated in Figure 2, showing plots of the response function R{t) against time t, both on 
a linear scale (left) and on a logarithmic scale over a larger range of times (right), for the 
same values of S. 





FIG. 2: Plot of the response function R{t) against time t for various values of 6. Left: linear scale. 
Right: logarithmic scale, over a larger range of times. 



VII. COMPLEX 'FISH BONE' SPECTRA IN THE NON-INTEGRABLE CASE 

Our earlier construction showing that the Markov matrix has an integrable spectrum 
with large degeneracies, is known to fail in the following circumstance: the parameter e is 
special (i.e., e = p/q with p + q < N — 1) and the zero-field rate reads w = ^ + 6 with 5 7^ 0. 

In this circumstance, the spectrum of the Markov matrix is indeed significantly different. 
Let us consider e = 1 and free boundary conditions for definiteness. 

The first occurrence of a non-integrable spectrum is at A^ = 3. The Markov matrix has 
the following 8 eigenvalues: 



A = 0, A = -l, A = -l 
X = -2-26± VI + 6 + 252, 
A = -2 - 2(5 ± ^5(25-1). 



2S, 



A = -2 - 2(5, 



(68) 



The last 2 eigenvalues are complex for < (5 < |, i.e., ^ < w < 1. 
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Similarly, for A^ = 4, the 16 eigenvalues of the Markov matrix read: 



A = 0, A = 

A = -1 - 25, 

A = -2- 

A = -3 

\ = -2 - 25 ± ^5{25 

A = -3 - 2(5 ± v^5(2^ 



1 (mult. 2), A = -2, 
A = -2-2(5 (mult. 2), 



A 



-3 - 2(5, 



2(5 ± Vl + (5 + 2(52, 
2(5 ± Vl + (5 + 2(52, 



(69) 



The last 4 eigenvalues are complex for < (5 < |, i.e., | < w < 1. The spectra however 
keep the inclusion property: the 16 eigenvalues (69) inclucie the 8 eigenvalues (68). 

As the system size A^ grows, more of the combinatorial degeneracies of the integrable case 
are lifted, and more eigenvalues become complex. These features are illustrated in Figure 3, 
showing the spectrum of the Markov matrix in the complex A plane, obtained by means 
of a numerical diagonalization, for e = 1 and A^ = 12, with free boundary conditions and 
for various values of 5. The spectra progressively become structured as 'fish bones' as |(5| 
increases. 
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FIG. 3: Plot of the spectra of the Markov matrices in the complex A plane, for A^ = 12 and 
various values of 5 (which are indicated to the right of the corresponding plot). Each spectrum is 
symmetric with respect to the real axis and consists of 2^2 = 4096 points. Spectra are translated 
vertically and shown in various colors for clarity. Left: (5 < 0. Right: (5 > 0. 



A quantitative measure of the extension of the fish bone spectra in the complex plane is 
provided by the mean squared imaginary part of the spectrum. 



K 



2 
a=l 



(70) 



This quantity is plotted in Fig. 4 against 5 = w — i^ for £ = 1, A^ = 10, and free and fixed 
boundary conditions. Data have been obtained by means of a numerical diagonalization of 
the Markov matrix for values of (5 on a grid with mesh 0.01. The integrable case [5 = 0) 
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is shown as a vertical blue line. The data suggest that k does not vanish in the vicinity of 
that point. In other words, the limits A^ — )■ oo and 5 — )■ do not commute. The quantity k 
rather exhibits a cusp (change of slope) at 5 = 0, whereas the observed dip is most certainly 
a finite- size effect. Finally, the mean squared imaginary part of the spectrum has a smooth 
maximum near 5 ~ 0.3, i.e., w ~ 0.8, irrespective of boundary conditions. 




FIG. 4: Plot of the mean squared imaginary part, «;, of the spectrum of the Markov matrix against 
8 = w — \ for £ = 1 and N = 10. Upper black curve: free boundary conditions. Lower red curve: 
fixed boundary conditions. Vertical blue line: integrable case ((5 = 0). 



A. Sensitivity: fish bones and false diagonalization 



In this section we study two related phenomena: for non-generic e, the rapidity with 
which the imaginary parts of the eigenvalues grow (as a function of 5), as exemplified in the 
"fish bone" graphs of Fig. 3, and another anomaly, mentioned in Remark 5, connected to 
numerical diagonalization in the case that a Jordan form is needed. What happens is that 
the computer spits out bona fide eigenvectors and non-degenerate eigenvalues, usually with 
small complex parts. The complex parts typically are much larger than the rounding error 
of the software. (This can also fool the investigator: if you get an imaginary part of order 
10^^ and if your accuracy is 10~^^, then you tend to believe the numerics. In our case we 
had analytic proof that the spectrum (with zero 5) consisted of real non-positive integers.) 

What we show in this section is how, for Jordan forms, their non-generic nature ex- 
presses itself in vastly exaggerated effects of rounding errors, including imaginary parts for 
real eigenvalues and rapid growth of the imaginary part when true perturbations move the 
spectrum off the real line. 

Consider the archetypal matrix requiring the Jordan form 



R = 



/O 1 ...\ 

10 ... 

1 

\0 000 00 0/ 



(71) 
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that is, the n-hj-n matrix, R, has non-zero entries only for R^i = Sk/+i. Of course this 
is already in Jordan form. Consider the effect of a single non-zero entry added somewhere 
to R. This could correspond to a rounding error (say of order 10"^^) or for the fish bone 
calculation this would be S. The additional term is xM, where x is a number and M is a 
matrix, all zeros, except for a single 1, whose location is described below. There are several 
cases. 

1. M's non-zero entry is on the diagonal. Then det (Al — R — xM) = A"^^(A — x). This 
implies that the shift in eigenvalue is just x, linear in the change. 

2. M's non-zero entry is anywhere above the diagonal. Then det (Al — i? — xM) = A", 
and there is no change in the eigenvalue. 

3. M's non-zero entry is on the j^^ sub-diagonal, that is Mki = 1 for k — i = j,n— 1> 
j > 0. Then det (Al — R — xM) = A"^-'^^ (A-'^^ — x). The non-zero eigenvalues are 
then {|a;|"'^'*--'^^''a;} where cu is one of the (j -|- 1)**^ roots of unity. 

The extreme situation occurs in Item 3, when x sits in the (n, 1) position, leading to n^^ 
roots of unity. For n > 2 these necessarily include complex values. For n = 2 the deviation 
can be real only. In any case it's clear how errors enormously larger than 10^^^ can occur. 
If there is a numerical error of the sort involving a Jordan block arising from degeneracy 
6 (which happens for > 4 spins) there can easily be complex roots of order 10~^. In an 
ordinary numerical calculation one would be hard put to call this a rounding error. 

In practice one does not start with a matrix of the form Eq. (71), but in the course of 
numerical diagonalization the many similarity transformations can move deviations of the 
pure Jordan form (whether intentional (6) or not (rounding error)), anywhere, including 
dangerous spots, like the (n, 1) position. Thus suppose the true Jordan form is achieved by 
the transformation w = V JV~^ with J the Jordan form and V the similarity transformation. 
Then a slightly modified w subjected to this same transformation will have small non-zero 
terms in "dangerous" locations, as described above. 

VIII. SUM RULE FOR THE MARKOV MATRIX 

In this section we provide an alternative investigation of the spectrum of the Markov 
matrix. We employ an explicit calculation of the average absolute eigenvalue 



2 



N 



A = -(A> = -^E^" = -^- P2) 



a=l 



This quantity can be evaluated by elementary means, thereby providing a sum rule for the 
spectrum. We have 

2iv 

^ = 4E^(C), (73) 



2^ 
c=i 



where C = {ai, . . . , a^} is an arbitrary spin configuration, and w(C) denotes the total exit 
rate from that configuration. 

Let us consider first free boundary conditions {hi = 0). For e generic, and for any 
configuration C = {ai, . . . , a^}, flipping the uppermost spin o"i always brings a contribution 
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w = ^ + S to oo{C), whereas the contribution coming from flipping any other spin (o"„ for 
n = 2, . . . ,N) is 1 if (Tn 7^ sign /i„ and if (T„ = sign /i„. We thus obtain 

A = - + <5. (74) 

For e = p/q special, the contribution to (^(C) coming from flipping o"„ is also w = ^ + 6 
whenever hn = 0. The number of configurations C such that /i„ = at depth n = k{p+q) + l 
reads 2^Pk, where 

P^ = -^ (k > 0) (75) 

is the probability for a random walker on the line, making integer steps +p and —q with 
equal probabilities, to be back to its starting point after k{p + q) steps. The result (74) is 
thus changed to 

A = ^ + A^S, (76) 

where for all A^ in the range k{p + q) + 1 < N < {k+ l){p + q), A]\f is constant and equal to 



Ak = Y,Pi- (77) 



1=0 

Let us now consider fixed boundary conditions {ctq = +1, hi = —1). For e generic, 

the contribution to uj{C) coming from fiipping any spin o"„ is 1 if (T„ 7^ sign/i„ and if 

(T„ = sign hn- We thus obtain 

A^ 
A = y. (78) 

For e = p/q special, the contribution to ti;(C) coming from fiipping (T„ is w = \ + S whenever 
hn = 0. The number of configurations C such that /i„ = at depth n = k{p + q) reads 
2^Qk, where 

^k{p+q)-l 

«' = ^k^ = JT7/' (* S 1) (79) 

is the probability for a random walker on the line, making integer steps +p and — g with 
equal probabilities, and starting from —q, to be at the origin after k{p + q) — 1 steps. The 
result (78) is thus changed to 

A^ 
A = y + Bm6, (80) 

where for all A^ in the range k{p + q) < N < {k + l){p + q) — 1, Bj^ is constant and equal to 

Bk = y^Qi = ^^{Ak-l). (81) 

The above results are summarized in the last column of Table III. In all the integrable 
situations where the spectrum is known explicitly, the corresponding simple values of A can 
be checked directly. The special results (76) and (80), with their non-trivial amplitudes Aj<4 
and i?7V; give some information on the non-integrable situations, namely the position of the 
center of mass of the corresponding fish bone spectra. 
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The generic results (74) and (78) can be recovered from the special ones (76) and (80) by 
taking the limit of an infinite period (p + g — )■ oo), so that Ajy = ^o = 1 and B^^ = Bq = 
for all system sizes N. Furthermore, if the integers p and q are interchanged, so that e is 
changed to its inverse, the P^ and the Af^ are left invariant, whereas the Qk and the Pat are 
multiplied by q/p. 

The case where e = 1, i.e., p = g = 1, is that of the usual symmetric walker. The 
probabilities 

P^ = Q^ = ^ (82) 

have a slow power-law decay, as Pk ~ (tt/c)^^/^. As a consequence, the amplitudes An and 
Bn grow as 

An^Bn- i—j . (83) 

The contributions proportional to S in (76) and (80) grow indefinitely, but sub-extensively, 
and therefore remain negligible with respect to the leading terms N/2. 

In all the other rational cases, i.e., e = p + q (irreducible) with p + q > 3, the walker is 
biased, and so the return probabihties decay exponentially, as Pk ^ Qk ^ a'^, with 

a = ^^^^ . (84) 

As a consequence, the amplitudes A^ and B^ saturate to finite limits: 

2p 

^ki -Oqo - 

r 

A;>0 



^oo = V Pk. B^ = ^(Aoo - 1). (85) 

^L^ p-\- q 



These limits turn out to be algebraic numbers. This property stems from the fact that the 
return probabilities P^ and Q^ are related to the Fuss-Catalan numbers [28]. Some variants 
of the above series have been met in several works in the physics literature [29-31]. 

To close with an example, for e = 2 we obtain A^o = 1 + S/y/b and B^o = 4/v^, whereas 
for e = 1/2 we obtain the same A^o, as expected, but B^o = 2/\/5. 



IX. BRINGING THE STOCHASTIC GENERATOR TO TRIANGULAR FORM 

In Remark 1 (in Sec. II B) we defined the matrix y (or y^) associated with a random walk 
on the edges of an A^-cube and presented its full spectral analysis. Its basis vectors, in Sec. 
V, acquired the interpretation of correlation functions, and we saw that in that basis the 
master equation has a triangular form. In this section we briefly derive that triangular form 
using the algebraic-recursion relation approach, noting as a consequence that its validity 
extends beyond the zero-temperature granular dynamics model. 

Let { V^} be the set of eigenvectors of y. Writing these vectors as the columns of a matrix 
V.a, we have X^^l/fc^^fo — ^aVka- The state numbering, and resultant labeling of V, follows 
the conventions for the matrix w. Then the matrix w = V~^wV is lower triangular. With 
this numbering, neither the eigenvalues of y nor those of w come in descending order. 
Proof: We first note the following recursion for the diagonalizing matrix for y. N.B. this uses 
the state-ordering scheme outlined in Sec. II A and is a direct consequence of the eigenvector 



construction. 
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"-Hi;: _;;:) t^) 



Moreover, it has the following properties: 



Recall too 



V;^ = Vm (87) 

V^' = VW/2^ (88) 



, Wiv - A A , 



We suppose that for A^, wn = Vj^'^wVn is lower triangular. This is the inductive hypoth- 
esis. It is true for A^ = 1 since in this case w coincides with y and is diagonalized by the 
transformation. 
Consider 

Performing the multiplications and making use of Eq. (88) yields 



Wn+I 



V^I^wn+iVn = ^ ( j_ 



VnWnVn 

+ 2KAK VnWnVn - 1 



Wjy 

-1 + 2-(^-iVA\/ wn-1 



(91) 



This expression is lower triangular. The exact form of A only enters in the lower left. 
This shows incidentally that all that is necessary for triangularity is the indicated recursion, 
recalling that there is also the condition A + A = 1. 

Remark 14: To obtain full agreement between the triangular matrix derived here and that 
obtained from the correlations, one must multiply each Va (as defined in Remark 1) by 
(— 1)'"^', with A the subset of {1,...,A^} associated with the index a. This is a unitary 
transformation of V. 



X. DISCUSSION 

Even with random shaking (positive "temperature"), the movement of grains has struc- 
tures and constraints not present in ordinary liquids. Without shaking, the zero-temperature 
dynamics has yet more restrictions and we have found the mathematical expression of this 
feature to be reflected in the matrix generator of the stochastic dynamics in the following 
way: for generic values of the shape parameter, e, the matrix has degenerate eigenvalues for 
which the number of eigenvectors is fewer than the degeneracy: it is not diagonalizable and 
the nearest one can get to a spectral expansion is the Jordan canonical form. 

In this article we have presented extensive information on the properties of the generator 
of the stochastic dynamics (which we call w). The mere fact that it obeys a particular 
recursion relation (our Eq. (7)) already fixes its eigenvalues and the multiplicity of the 
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associated invariant spaces [invariant space is what replaces the notion of eigenvector span 
when they do not span). Moreover, the recursion is more general than the granular model 
in the sense that other matrices also satisfy the recursion (they need not even be stochastic) 
and have the same spectrum. Nevertheless, this same spectrum was derived in a way that 
made direct use of the model's features by studying correlation functions among the spins. 
This approach made use of a fundamental asymmetry in the model: a given spin is affected 
only by those above it. The triangular structure of the differential equations obtained in 
this way is identical that obtained from a particular similarity transform applied to the 
stochastic generator. We also show that the recursion property alone is enough to guarantee 
that this transform puts w into (lower) triangular form. 

The particular similarity transform casts light on one of the lovely features of the model, 
namely the appearance of combinatorial coefficients for the multiplicity of the invariant 
spaces. In particular, the eigenvalues of w are {0, —1, —2, . . . , — A^} with A^ the number of 
spins in the vertical column that models the granular material. The multiplicity of eigenvalue 
—k is C^ (we take S = in this discussion). Now the invariant spaces are not character- 
ized by the number of spins in one direction or another, which would have been a natural 
way to obtain combinatorial coefficients of the sort indicated. However, the triangulating 
transformation does count spins in the following sense. It is the diagonalizing matrix of "y," 
described in Remark 1. That matrix also describes A^-spin stochastic dynamics, but for it, 
all single spin ffips are allowed, in contrast to the granular model with its highly structured 
rules for transitions. The eigenvectors of y of eigenvalue —k (and it has a full complement, 
unlike w) are obtained by looking at the subsets of size k, so indeed there will be C^ of 
them for each k. These subset-based eigenvectors are what triangularize w. 

The relation of the eigenvalues and the combinatorial coefficients is also clear in the 
development of the differential equations for the correlation functions (Sec. V). When each 
spin can ffip independently, the combinatorial coefficients, as well as the full spectrum, can 
be obtained by looking at subsets of the spins of particular cardinality — equivalent to the 
diagonalization of the generator of the random walk on the hypercube (what we called y). 
For the model discussed here we do have dependence on the motion of other spins, but 
it is asymmetric, with a given spin's dynamics depending only on those above it. This 
is sufficient to imply that the differential equations for the correlations have a triangular 
form, and that both the eigenvalues and the cardinality of the invariant spaces is unchanged 
from the all-independent model. Finally, in a satisfying conffuence of the two approaches, 
all coefficients in the triangular form obtained in this way coincide with those of the pure 
matrix approach. 

Our mathematical developments have not completely answered, "why Jordan?" nor 
the yet more compelling question, "why maximally Jordan?" Certainly this property is a 
reffection of the many constraints in the system. These constraints are in their own way 
"maximal," in the sense that for every potential spin ffip at level n, all spins above n are 
queried. It is thus plausible that this maximal constraint on conffgurations induces the 
cascading structure in the "maximally Jordan" invariant spaces. In the time domain, a 
consequence of the Jordan form is that the relaxation process is not direct, and we have 
explored the way in which the Jordan form induces non-monotonic behavior, despite the 
presence of entirely real spectrum (for non-generic e and nonzero 6 complex eigenvalues 
obviously allow non-monotonicity). In Sec. IIIB we gave a mathematical criterion for the 
maximal Jordan property and from the correlation function approach the interdependence 
and triangular form of the matrix suggests this property. But it must be cautioned that 
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the recursion alone, or triangularity alone does not compel the maximal Jordan property 
and we have produced examples (not presented in this article) in which various patterns 
of the matrix A^v (diagonal and with O's I's and 1/2's on its diagonal), when used in the 
recursion, do not lead to the maximal Jordan property. Nevertheless, at the heuristic level 
our picture is that a cascade of processes takes place, with earlier transitions needed before 
later ones can happen. This is different from what happens in a (say) spin glass where 
individual configurations must do unlikely things in order to move to lower energy states. In 
the grain model it is the probability distribution that behaves non-monotonically, which is 
not the case for the spin glass. However, the "cascading" is in a way less dramatic than what 
happens in the spin glass. There, arriving at a lower energy state can take exponentially long; 
for us, the delay only slows things down in a polynomial way, with the dominant (and not 
particularly slow) exponential ultimately manifesting itself. Taking heed of this example, we 
note that "hindered," "complex dynamics" and "Jordan forms" are independent concepts. 
In particular, most models exhibiting glassy dynamics and aging have explicit reference to a 
static Hamiltonian. The resulting stochastic dynamics therefore obeys detailed balance. The 
associated Markov matrix can be brought to a symmetric form. It is therefore diagonalizable, 
and its spectrum is real. These features are best exemplified by the so-called "kinetically 
constrained" models [32], such as the Frederickson- Andersen spin model [33], where entropic 
barriers in the configuration space place restrictions on allowed transitions, but, as indicated, 
do not require the Jordan form. As a corollary, in these models the time dependence of the 
probability distribution is always monotonic in its approach to the stationary (equilibrium) 
state, a property we have explicitly shown not to hold for our model. 

The principal results of the present work are concerned with the short-time relaxation 
of the model, dealing as they do with eigenvalues and with the huge — for large magnitude 
eigenvalues — invariant spaces of the Markov matrix. This does not deal with many other 
physical issues, many of which are concerned with the structure of the stationary state. Thus 
many of the short time properties we have found hold irrespective of boundary conditions 
(see Table III), whereas the essential physical attributes of the model are sensitive to this 
feature [9-12]. For free boundary conditions (the situation implicitly under consideration 
until Sec. VF) the stationary state is always a non-trivial fluctuating non-equilibrium steady 
state. For fixed boundary conditions (e.g., (Tq = +1), the nature of the stationary state 
depends on parameters. For e generic (i.e., essentially irrational), the system is driven to a 
unique ground state. The ground-state configuration is absorbing for the zero-temperature 
dynamics. For an infinite system, it is a quasiperiodic sequence, which admits a geometrical 
construction. For e special (i.e., essentially rational), so that there are internal points of 
zero field, the system reaches a fluctuating stationary state characterized by an anomalous 
scaling of fluctuations. 
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Appendix A: Conjugation transformation 



In this Appendix we perform — at the matrix level — the transformation carried out on the 
determinant in Eqs. (10)-(14). It is a systematic procedure, different from that appearing 
in Sec. IX, for obtaining a triangular form of the matrix. 

As above, 1^ and 0^ are the 2^x 2^ unit and zero matrices, respectively. Let a^ = 

-__;^T, ; _-i Now fix A^ (the number of spins). We next 



0;. li 



, with inverse a^^ 



0. 



\(N) 



define A^ (or Ak, with the N suppressed). For this N, concatenate 2 a^'s along diagonal 

blocks so as to make a 2^x2^ matrix. This concatenated object is A^ . Thus for k = 1, 
ai is a 2x2 matrix, so that Ai has 2^^^ such diagonal blocks: 



Ai 



( «i Oi 
Oi ax 

VOi Oi 



Oi\ 



Ol / 



(AI) 



For given A, the largest "A^" is A^^x (which is composed of two copies of a^). The inverse 

a f3 
7 6 
7 —a — 7 + /3 + 5 



of Ak is constructed in the same way from the inverses of a^. If Q 



then 



a 



. For our matrices, either w^ or w^- — AI 



N, 



aQ a ^ . 

7 -7 + 

the upper right block of P is zero, so that a succession of A transformations can bring w 
to lower triangular form. The procedure for doing this is to work your way down: First 
transform wn by Ajsf-i then by A7V-2; and so forth, with the final transformation using Ai. 
What then sits on the diagonal are the eigenvalues of w. 

This technique works whether or not a Jordan form is required. Thus for the fully 
connected graph the recursion that yields the continuous time generator uses A = (1/2)1 at 
each stage of the construction. As for w, the associated matrix is disassembled by the A's, 
but since there are 2^ eigenvectors, no Jordan form is needed. 



Appendix B: Quasisymmetries 



1. The operator F 



Recall that the (diagonal) matrix in the lower- left block of w is A (see Eq. (7)). Although 
it might seem that the and 1 count should be the same for both it and A, there is 
nevertheless an important difference. The role of source and target is reversed and therefore 
the sequence is reversed. This can be formally stated by defining a family of matrices F/v 
to be the left-right reversal of the 2^x2^ identity (Iat). Thus for example 



/O 1\ 
10 
10 

Vi 0/ 



(Bl) 
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Note that F^ = 1. With this notation 

An = FnA^F^, ViV. 



(B2) 



2. Intertwining operators 



Let 



which is a mapping from 



N 



[1,1] 



Lat 



(B3) 



d2^+i 



— )> 



P, 



and (e.g.) for A^ = 2 takes the form 

/I 1 0\ 

10 10 

10 10 

VO 1 1/ 



(B4) 



Ahernatively, P^ 



LJV 



LAT 



Interest in this operator arises from the identity 



2A 



N 



Pn ^Af+l Pn 



(B5) 



where An can be either wn or its transpose (^at. This operator has an interpretation as a 
projection, mapping pairs of spin states onto a single spin state. It satisfies the identities 



PnPn 



21 



N 



Pn^n 



1 1 
1 1 



lat 



(B6) 



The "doubling" of Sec. Ill is multiplication by P^. Also note (which is what lies behind the 
demonstrations of that Section) that 



w 






Af+1-' N 



W 



N 



W 



N 



(B7) 



since the A's are eliminated, as in Eq. (24). 

The operators P and F account for the symmetries observed in the form of the eigenvec- 
tors. 
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